home *** CD-ROM | disk | FTP | other *** search
/ Aminet 21 / Aminet 21 (1997)(GTI - Schatztruhe)[!][Oct 1997].iso / Aminet / gfx / x11 / xmountains.lha / src / random.c < prev    next >
C/C++ Source or Header  |  2003-01-20  |  4KB  |  184 lines
  1. /*
  2.  *    C version of Marsaglia's UNI random number generator
  3.  *    More or less transliterated from the Fortran -- with 1 bug fix
  4.  *    Hence horrible style
  5.  *
  6.  *    Features:
  7.  *        ANSI C
  8.  *        not callable from Fortran (yet)
  9.  */
  10.  
  11.  
  12.  
  13. char uni_id[] = "$Id: random.c,v 1.6 1995/01/20 15:13:06 spb Exp spb $" ;
  14. /*
  15.  *    Global variables for rstart & uni
  16.  */
  17.  
  18. #define PARANOID
  19. /* need types and time */
  20. #include <X11/Xos.h> 
  21. /* #include <sys/types.h>
  22.  * #include <sys/time.h>
  23.  */
  24. #include <stdio.h>
  25. #include <math.h>
  26. typedef struct
  27. {
  28.   float u[98];
  29.   float c;
  30.   float cd;
  31.   float cm;
  32.   int ui;
  33.   int uj;
  34. }
  35. Uni_save;
  36.  
  37. Uni_save uni_data;
  38.  
  39. float uni()
  40. {
  41.     float luni;            /* local variable for uni */
  42.  
  43.     luni = uni_data.u[uni_data.ui] - uni_data.u[uni_data.uj];
  44.     if (luni < 0.0)
  45.         luni += 1.0;
  46.     uni_data.u[uni_data.ui] = luni;
  47.     if (--uni_data.ui == 0)
  48.         uni_data.ui = 97;
  49.     if (--uni_data.uj == 0)
  50.         uni_data.uj = 97;
  51.     if ((uni_data.c -= uni_data.cd) < 0.0)
  52.         uni_data.c += uni_data.cm;
  53.     if ((luni -= uni_data.c) < 0.0)
  54.         luni += 1.0;
  55.     return ((float) luni);
  56. }
  57.  
  58. void rstart(i,j,k,l)
  59. int i;
  60. int j;
  61. int k;
  62. int l;
  63. {
  64.     int ii, jj, m;
  65.     float s, t;
  66.  
  67.     for (ii = 1; ii <= 97; ii++) {
  68.         s = 0.0;
  69.         t = 0.5;
  70.         for (jj = 1; jj <= 24; jj++) {
  71.             m = ((i*j % 179) * k) % 179;
  72.             i = j;
  73.             j = k;
  74.             k = m;
  75.             l = (53*l+1) % 169;
  76.             if (l*m % 64 >= 32)
  77.                 s += t;
  78.             t *= 0.5;
  79.         }
  80.         uni_data.u[ii] = s;
  81.     }
  82.     uni_data.c  = 362436.0   / 16777216.0;
  83.     uni_data.cd = 7654321.0  / 16777216.0;
  84.     uni_data.cm = 16777213.0 / 16777216.0;
  85.     uni_data.ui = 97; /*  There is a bug in the original Fortran version */
  86.     uni_data.uj = 33; /*  of UNI -- i and j should be SAVEd in UNI()     */
  87. }
  88.  
  89.  
  90. /* ~seed_uni: this takes a single integer in the range
  91.  *        0 <= ijkl <= 900 000 000
  92.  *    and produces the four smaller integers needed for rstart. It is
  93.  *    based on the ideas contained in the RMARIN subroutine in
  94.  *        F. James, "A Review of Pseudorandom Number Generators",
  95.  *            Comp. Phys. Commun. Oct 1990, p.340
  96.  *    To reduce the modifications to the existing code, seed_uni now
  97.  *    takes the role of a preprocessor for rstart.
  98.  *
  99.  */
  100.  
  101. void seed_uni(ijkl)
  102. int ijkl;
  103. {
  104.     int i, j, k, l, ij, kl;
  105.  
  106.     if( ijkl == 0 )
  107.     {
  108.       ijkl = time((time_t *) 0);
  109.       ijkl %= 900000000;
  110.     }
  111.     /* check ijkl is within range */
  112.     if( (ijkl < 0) || (ijkl > 900000000) )
  113.         {
  114.         fprintf(stderr,"seed_uni: ijkl = %d -- out of range\n\n", ijkl);
  115.         exit(3);
  116.                    }
  117.  
  118.  
  119.     /* decompose the long integer into the the equivalent four
  120.      * integers for rstart. This should be a 1-1 mapping
  121.       *    ijkl <--> (i, j, k, l)
  122.      * though not quite all of the possible sets of (i, j, k, l)
  123.      * can be produced.
  124.      */
  125.  
  126.     ij = ijkl/30082;
  127.     kl = ijkl - (30082 * ij);
  128.  
  129.     i = ((ij/177) % 177) + 2;
  130.     j = (ij % 177) + 2;
  131.     k = ((kl/169) % 178) + 1;
  132.     l = kl % 169;
  133.  
  134. #ifdef PARANOID
  135.     if( (i <= 0) || (i > 178) )
  136.         {
  137.         fprintf(stderr,"seed_uni: i = %d -- out of range\n\n", i);
  138.         exit(3);
  139.                    }
  140.  
  141.     if( (j <= 0) || (j > 178) )
  142.         {
  143.         fprintf(stderr,"seed_uni: j = %d -- out of range\n\n", j);
  144.         exit(3);
  145.                    }
  146.  
  147.     if( (k <= 0) || (k > 178) )
  148.         {
  149.         fprintf(stderr,"seed_uni: k = %d -- out of range\n\n", k);
  150.         exit(3);
  151.                    }
  152.  
  153.     if( (l < 0) || (l > 168) )
  154.         {
  155.         fprintf(stderr,"seed_uni: l = %d -- out of range\n\n", l);
  156.         exit(3);
  157.                    }
  158.  
  159.     if (i == 1 && j == 1 && k == 1)
  160.         {
  161.                 fprintf(stderr,"seed_uni: 1 1 1 not allowed for 1st 3 seeds\n\n");
  162.         exit(4);
  163.                 }
  164. #endif
  165.  
  166.  
  167.         rstart(i, j, k, l);
  168.  
  169. }
  170.  
  171. float gaussian()
  172. {
  173.         double pi = 3.1415926536, two = 2.0, zero = 0.0;
  174.     double ran1, ran2;
  175.  
  176.     do {
  177.         ran1 = (double) uni();
  178.     } while (ran1 == zero);
  179.  
  180.     ran2 = (double) uni();
  181.     return (float) ( sqrt(-two * log(ran1)) * cos(two * pi * ran2) );
  182. }
  183.  
  184.